Binomial test and binomial confidence intervals with python

Author: Leonardo Espin

Date: 8/26/2019

In [1]:
import numpy as np
from math import factorial #for binomial coefficient
from scipy.stats import norm #for normal approximation of distribution of binomial proportions
from scipy.stats import binom #for binomial distribution
In [2]:
def binomial(i,n,p):
    '''
    Calculation of the binomial distribution using the formula 
    (not good for large values of n):
    '''
    return (factorial(n)/(factorial(i)*factorial(n-i)))*(p**i)*((1-p)**(n-i))

  • Example: A four-sided (tetrahedral) die is tossed 1000 times, and 290 fours are observed. Is there evidence to conclude that the die is biased, that is, say, that more fours than expected are observed?

  • We count events as rare as, and rarer than the observed proportion, meaning 290 tosses or more

In [3]:
#direct calculations of probability of events as rare as, and rarer than 
#the observed proportion, i.e 290 or higher events
print('p-value for observed proportion: {:.4f}'.format(
    sum([binomial(i,1000,0.25) for i in range(290,1001)])))

p-value for observed proportion: 0.0022

Therefore there is significant evidence (p < 1%) that the null hypothesis (dice is fair) is false.

We could also have calculated this by using the binomial survival function (1-cdf, the cumulative density function):

In [4]:
print(sum([binomial(i,1000,0.25) for i in range(290,1001)]))
print(binom.sf(289,1000,0.25)) #for some reason, the survival function starts summing from i+1 not i?

0.0022296751567737897
0.0022296751567762976
  • Since the sample is large (n=1000) it is advisable to use the normal approximation via the z-statistic (we are calculating the binomial coefficient for large values)
    • For example, the 0.05 p-value for a 0.25 proportion in a sample of a thousand is:
In [5]:
#using the inverse survival function of the binomial distribution
binom.isf(0.05,1000,0.25)
Out[5]:
273.0

now let's calculate the p-values for a proportion of 0.273:

In [6]:
print(sum([binomial(i,1000,0.25) for i in range(273,1001)]))
print(binom.sf(272,1000,0.25)) #for some reason, the survival function starts summing from i+1 not i?

0.0511946713027778
0.05119467130280299

now with the normal approximation:

In [7]:
norm.sf((0.2725-0.25)/np.sqrt(0.25*0.75/1000))#notice that we requiere more accuracy, 0.2725 instead of 0.273
Out[7]:
0.0501741232311452

or equivalently:

In [8]:
norm.sf(0.2725,loc=0.25,scale=np.sqrt(0.25*0.75/1000))#p-value, the mean and standard deviation respectively
Out[8]:
0.0501741232311452

This interval certainly contains the population mean. However, notice that it is contained in the bootstrapped confidence interval above, and we don't know what is the probability attached to this interval of it containing the mean.

In [9]:
norm.isf(0.05,loc=0.25,scale=np.sqrt(0.25*0.75/1000))
Out[9]:
0.2725230858818877

  • Example: Among patients with lung cancer, usually 90% or more die within three years. As a result of new forms of treatment, it is felt that this rate has been reduced. In a recent study of n = 150 lung cancer patients, y = 128 died within three years. Is there sufficient evidence at the α = 0.05 level, say, to conclude that the death rate due to lung cancer has been reduced?

  • Note that here counting events as rare as, and rarer than the observed proportion means events of 128 deaths or less

In [10]:
#counting events of 128 deaths OR LESS
print('p-value for observed proportion: {:.4f}'.format(
    sum([binomial(i,150,0.9) for i in range(128)])))

#using the binomial cumulative density function
print('p-value for observed proportion: {:.4f}'.format(
    binom.cdf(127,150,0.9)))#again, the survival function starts summing from i+1 not i?

#using the normal approximation
print('p-value for observed proportion: {:.4f}'.format(
    norm.cdf((128/150 -0.9)/np.sqrt(0.9*(1-0.9)/150))))

p-value for observed proportion: 0.0256
p-value for observed proportion: 0.0256
p-value for observed proportion: 0.0284

Therefore we can reject $H_0$ = "the rate is still 0.9" at the α = 0.05 level. Note, however, that we can't reject $H_0$ with 1% confidence.

Indeed, the event of a proportion of 128 out of 150 is barely within the 0.05 confidence band:

In [11]:
#using the percent point function (inverse of cdf — percentiles)
binom.ppf(0.05,150,0.9)
Out[11]:
129.0

Calculating a 95%, symmetric confidence interval with the normal approximation:

In [12]:
#Inverse survival function (inverse of sf).
print('upper-bound proportion: {:.4f}'.format(
    norm.isf(0.025,loc=0.9,scale=np.sqrt(0.9*(1-0.9)/150))))#loc is the mean, scale the st. dev.

#Percent point function (inverse of cdf — percentiles)
print('lower-bound proportion: {:.4f}'.format(
    norm.ppf(0.025,loc=0.9,scale=np.sqrt(0.9*(1-0.9)/150))))

print('observed proportion: {:.4f}'.format(128/150))

upper-bound proportion: 0.9480
lower-bound proportion: 0.8520
observed proportion: 0.8533

So, the observed proportion is inside the 95% symmetric confidence interval and we can't reject $H_0$